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We study analytically and numerically the corrections to scaling in turbulence 
which arise due to the finite ratio of the outer scale L of turbulence to the viscous 
scale 7], i.e., they are due to finite size effects as anisotropic forcing or boundary 
conditions at large scales. We find that the deviations 5( m from the classical Kol- 
mogorov scaling ( m = m/3 of the velocity moments (\u(k)\ m ) oc k~^ m decrease 
like 5( m (Re) = CmRe^ 10 . Our numerics employ a reduced wave vector set ap- 
proximation for which the small scale structures are not fully resolved. Within this 
approximation we do not find Re independent anomalous scaling within the inertial 
subrange. If anomalous scaling in the inertial subrange can be verified in the large 
Re limit, this supports the suggestion that small scale structures should be respon- 
sible, originating from viscosity either in the bulk (vortex tubes or sheets) or from 
the boundary layers (plumes or swirls). 

PACS: 47.27.-i, 47.27.Jv, 47.27.Eq, 47.27.Nz 
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A major question in the theory of turbulent flows is whether there exists an 
asymptotic scaling state at high Re, and whether in this state the values of the 
scaling exponents conform with the classical predictions of the Kolmogorov the- 
ory [[I]]. There is some experimental evidence that a scaling state exists, and that 
the scaling exponents deviate from the Kolmogorov prediction (anomalous scaling). 
Theoretically there is still no proof for anomalous scaling based on the Navier-Stokes 
equations. Recently it was suggested that the reason for deviations from the Kol- 
mogorov theory is related to the creation of small scale structures like vorticity 
sheets and tubes in which the local geometry is not 3-dimensional and isotropic M. 

One could hope that numerical simulations could be used to decide whether 
classical or anomalous scaling should be expected in high Re flows. If anomalous 
scaling were found, simulations could distinguish the responsible physical mech- 
anism. Unfortunately, in direct simulations of the Navier-Stokes equations with 
current computers one can achieve only up to Taylor-Reynolds numbers of about 
Re\ = 200 ||. The range of scales for which power law behavior is observed at such 
values of Re\ is too small to distinguish between classical and anomalous scaling. 
In addition, at low to moderate values of Re one can have corrections to scaling due 
to finite size effects, and it is hard to take these into account if their expected Re 
dependence is not known. It is this latter issue which is the focus of this letter. We 
examine these corrections to scaling analytically and numerically, and show how to 
take them into account in any future experiment or numerical simulation. 

A high Re flow can be constructed in a numerical simulation by using a reduced 
wave vector set approximation (Fourier- Weierstrass decomposition), that was intro- 
duced and studied extensively by two of us recently. For a detailed description and 
references of previous work we refer to refs. 0, [5], || . The main idea is to start with 
a regular Fourier decomposition of the velocity field u(x, t) in terms of plane waves 
exp (ik ■ x), but admit only a geometrically scaling subset K = U[K t of wave vectors 
in the Fourier sum, 

u(x, t) — u(k, t) exp (ik ■ x). (1) 

k&K 

The subset of wave vectors is chosen such that Kq = {k^\n = 1, . . . , N} and 
Ki = {2 l k^\ n = 1, . . . , N}, 1 = 1,..., l max . The basic set K is chosen to contain 
wave vectors of different lengths that interact dynamically to a good degree, and 
lmax is chosen large enough to guarantee that the amplitudes u(k^ max \ t) of the 
smallest scales are practically zero. This of course depends on the viscosity v and 
the Reynolds number Re. The flow is driven by a deterministic, non stochastic 
driving f(k,t) with k in Kq. In fact, only the seven smallest wave vectors in K 
were driven. Obviously, this driving is not isotropic. The fc-range of externally 
forced amplitudes is denoted henceforth as stirring subrange. 

The simulations that we refer to here used a value of N = 80 and scanned a 
range of Re from 10 4 to 10 7 |J. It should be stressed that in this scheme the density 
of wave vectors per k interval decreases like 1/k, whereas it increases like k 2 in full 
grid simulations. For this reason small scale structures (as vortex tubes or sheets) 
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cannot be resolved very accurately. But on the other hand many more scales than 
in full simulations can be examined. For the largest Re our simulations comprise 
more than three orders of magnitude between the wave number with maximum 
dissipation and the outer scale. The scaling exponents ( m can be determined with 
high accuracy. 

The scaling exponents were measured in this numerical simulation by computing 
the m -order spectra (\u(k)\ m ), and fitting the three parameter function 

(\u(k)\ m ) = c m k-<-exp[-k/k i r ) }. (2) 
It was found [!J that a good approximation could be obtained globally with 

kt ] = 2kf/m, kf = (lS.Sr?)- 1 , (3) 

with 7] being the Kolmogorov length scale. Eq. (3) means that for higher order 
velocity moments the crossover between the inertial range and the viscous range 
takes place at smaller value of k, i.e., the inertial range is less extended. As we want 
to examine the effects of the large scales, i.e., of the stirring subrange, we eliminate 
viscous effects by fitting (2) with fixed k^ (according to (3)) only in the interval 
[0, kj]. The resulting scaling exponents are (very slightly) corrected to guarantee 
(3 = 1, see m H for the procedure. (3 = 1 is required by Kolmogorov's structure 
equation J7J, which strictly holds only in an isotropic situation. But as (3 = 1 is 
frequently enforced in the analysis of experimental data, we do it here as well. We 
represent the numerical results for the scaling exponents in terms of the deviations 
from the classical predictions ( m (Re) — m/3. The results, which pertain now to 
both inertial range and stirring range, are shown in fig. 1 for four different Re. The 
deviations are decreasing with increasing Re. 

Qualitatively this feature can be understood from the comparison with the local 
5(m(k) for each Re, see ref. ||. The S( m (k) are calculated by locally employing 
the fit (2) to the spectra with k^ fixed. The results are shown in fig. 2a for 
Re = 1.05 • 10 4 and in fig. 2b for Re = 1.4 • 10 7 . Clearly, with increasing Re the 
inertial range extension exceeds more and more that of the stirring range, and thus 
the effect of the nonuniversal forcing diminishes. 

To get quantitative information, we fitted the power law 

6( m (Re) = c m Re~^ (4) 

to the 5( m (Re) determined from the spectra. The results of the fit (4) for c m and 
j3 m are given in table 1. We find that (3 m is close to 3/10 for all m. 

A good way to understand this behavior is to consider the Euler part of the 
equations of motion using the representation by the Clebsch variables. In these 
variables the velocity field is written in terms of the two scalar functions X(x,t) and 
H(x,t): 

u(x,t) = XV /j, - V0, (5) 
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where the potential 4>(x, t) is determined from the incompressibility condition. Using 
these variables, the equation of motion 

d t u(x, t) + u ■ Vu = — Vp(a3, t) (6) 

reads 

8T~L 51~~L 

d t X{x,t) = d t fi(x,t) = -— — — , (7) 

dfi(x,t) dX(x,t) 

where the Hamiltonian 7i is defined as Ji = J dx\u(x, t) | 2 /2. This formulation allows 
to introduce normal canonical variables a(x,t) and a*(x,t), cf. ||, 

V2a(x,t) = X(x,t)+ifi(x,t). (8) 

These canonical variables can be used to expose integrals of motion. In particular, 
the "occupation number" N and the "momentum" P of quasi particles, 



N 



J \a k \ 2 dk, P = Jk\a k \ 2 dk, (9) 



are such integrals of motion [^, [IU|. These integrals of motion could also be ex- 
pressed as functions of the velocity field u(x,t), but are highly nonlocal in that 
representation. The main point is that the integral P can be non zero only in an 
anisotropic system, and therefore can lead us to scaling laws that capture the rate of 
decay of anisotropy as a function of k. To proceed in this direction we first find by 
dimensional analysis that in an isotropic system the k-component of N, N k , depends 
on k and the mean energy flux per unit time per unit mass, e, according to 

N k 5{k - k') = (a k a* kl ) = Ce 1/3 k- 13/3 5{k - k') (10) 

with C being a dimensionless constant. Of course, this result can be translated back 
to the standard Kolmogorov result for the energy E k = Ce 2 ^ 3 k~ 5 ^ 3 . Next, we realize 
that in an anisotropic system this relation can change, since N k can depend now on 
a dimensionless ratio of the "momentum flux" over the energy flux e. The energy 
flux is defined by the equation 



d t H k + -^ = 0, (11) 



dk 



which in a stationary state requires e k to be fc-independent, e k = e. The momentum 
flux 7Tfc, likewise, is defined by 

- d t P k = kk 2 d t N k (t) = (12) 

In a stationary system 7r k is also /c-independent, and from eqs. (0) and (0) we see 
that the dimension of tt can be written as [n] = [e 2 / 3 ^ 1 / 3 ] = length/ (time) 2 . In an 
anisotropic system we seek a new solution for Nk, 

N k = Ce^k-^fHk) (13) 
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where the dimensionless parameter is proportional to the ratio of the two fluxes tt 
and e. Since the two fluxes e and n have different dimensionalities, the parameter 
must involve k and e to some powers. The unique combination that is proportional 
to Ti is 

7T 

£ fc = e 2/3fcl/3 ' 

The function /(£&) is chosen such that for n = Eq. (0) regains its Kolmogorov 
form, i.e. /(0) = 1. Asuming that for small 7r the correction to the spectrum is 



proportional to tt, for small £ we choose /(£) = 1 + 0(£). Notice that in C p~3|) and 
(|TJJ) we assumed that the anisotropic correction to N k is analytic in tt. This is of 
course a crucial assumption, but it can be justified by an explicit calculation, and it 
was shown to be correct to all orders in perturbation theory, see ||. Denoting the 
anisotropic correction as 5N k , we conclude that for small anisotropy, 1, 

5N k /N k oc fc~ 1/3 . (15) 

Obviously, since TT k is odd in k, cf. (|D|) or (|T2D, so is also 5N k . In contrast, the 
double correlation function of the velocity field F(k), 

F{k)A(k- k') =Tr(u{k)u*{k')) (16) 

is even in k. (Note that in accordance with our discrete simulations we define F(k) 
with a Kronecker A.) Now, since F(k) is even in k, the lowest order correction due 
to anisotropy must be the second order. We can conclude therefore that 

5F{k)/F(k) oc k~ 2/s . (17) 

To bring this result to a form that can be tested against the simulation, we multiply 
F(k) by k (as numerically we of course use discrete wave vectors) and rewrite our 
result in the form 

F(k) = (\u(k)\ 2 ) = c 2 AT 2 / 3 (l + a m (kL)- 2 ^) exp \-k/kf\. (18) 

(L is the outer length scale.) The same conclusion is obtained by considering the 
effect of external shear |IT| . Since it introduces another time scale r(ki), spectral 
corrections should be oc r(k) / t^l) oc k~ 2 ^ 3 . On scales on which the shear dominates 
(k <^ L^ 1 ) the second term in ([18]) will be even the leading contribution and thus 
will change the dominating /c-power to fc _4//3 . For a similar recent result concerning 
the fc _4//3 -scaling law for small k in an anisotropic flow see [I2|. The crossover 



from k 4//3 to k 2//3 can well be observed in experiments [12|. Here we restrict 
ourselves to k > L^ 1 . So the second term in (JIB|) is an anisotropy correction, which 



will become smaller for increasing k due to isotropization by eddy decay. Yet it 
will render the scaling exponent k-dependent. Denoting the local scaling exponent 
dln(\u(k)\ 2 ) I dink as ^(k), we find from ( Jl8| ) 

C2W = I + w + 3[i+a,-W 3 r (19) 
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As a function of k this expression has a minimum which is rather flat, and can 
lead to an apparent exponent Q° ■ We estimate the value of this exponent by the 
minimum of (|19|) which approximately is 

d „> ' + -( » ] 2/5 W 3/ 5 . (20) 




Here, we used the fact that l<a 2 l (kL) 2 ^ 3 near the minimum, justified a posteriori, 
as we can predict that the deviation of ^ app ' } from 2/3 

goes down when Re increases, since kfh oc Re 3 / 4 0. We obtain 

5d app) oc Re- 3/1 °. (21) 

Exactly the same calculation can be repeated for all other scaling exponents ( m , and 
all of them approach their Kolmogorov 41 values with the same scaling oc Re~ 3 ^ 10 . 
The result of the calculation is 

ti PP) = J + *M<*n)™ (^) 2/5 K! 3/5 , (22) 

SCt Pp) = c m Re- 3 ' w (23) 

with 

^ = s ^(^)(™)2/5|^|3/5_ (24) 

c 2 a 2 2 a 2 

From eq. (|4]) and table 1 we can determine the m-dependence of a m . Assuming 
a power law dependence, we estimate a m oc 77i 3 - 8±0 - 8 ; which is a rather strong in- 
dependence. This finding for a m means that higher order moments are more affected 
by corrections to scaling than lower order ones. Together with the fact that their 
power law behavior is cut off at lower k vectors, cf. eq. (3), it means that it is harder 
to observe clean power laws in higher order moments for an appreciate range of 
scales. 

In summary, the results of the numerical simulations of the Navier-Stokes equa- 
tion with a reduced wave vector set approximation are in very close correspondence 
with the predictions of a theory in which the deviations from classical scaling is 
solely due to corrections to scaling. It is tempting to conjecture that if a full treat- 
ment of the Navier-Stokes dynamics results in genuine anomalous scaling that does 
not disappear in the limit of high Re, it has to do with the small scale structures 
that are underestimated in our simulations as well as in the analytical argument 
that refered to the Euler equations and neglected viscosity and boundary layers. 
Such a conjecture would be in accord with other theoretical expectations which link 
anomalous scaling to the creation of vortex structures like tubes and sheets. 
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Tables 



m 


2 


4 


6 


8 


10 




0.24 


0.29 


0.28 


0.31 


0.30 


Cm 


+0.11 


-0.32 


-1.30 


-4.03 


-7.02 



Table 1 

Fitparameters to eq. (4). If the m-dependence of c m is expressed by a power law, we 
will get c m oc m 2 - 7±0 - 5 . The data are compatible with c m oc (m — 3) 3±0 5 , indicating 
agreement with the Kolmogorov structure equation, which forbids corrections if 
m = 3. Here also the change of sign of the c m is immediately grasped. 
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Figure 1: Results from our reduced wave vector set approximation for \S( m (Re)\ 
for m=2,4,6,8,10, bottom to top. The straight lines correspond to the scaling law 
5( m (Re) oc Re' 3 / 10 as predicted by fl23|). 



Figure 2: Scale resolved intermittency corrections — S( m (k) for m = 2,4,6,8,10, 
bottom to top. In (a) we have Re = 1.05 • 10 4 , in (b) it is Re = 1.4 ■ 10 7 . The fit 
range is [k/y/lO, k^/lO] for all k. For details see ref. 0. 0. 
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